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ABSTRACT 

Rotational shear layers at the boundary between radiative and convective 
zones, tachoclines, play a key role in the process of magnetic held generation in 
solar-like stars. We present two sets of global simulations of rotating turbulent 
convection and dynamo. The hrst set considers a stellar convective envelope 
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only; the second one, aiming at the formation of a tachocline, considers also 
the upper part of the radiative zone. Our results indicate that the resulting 
mean-flows and dynamo properties like the growth rate, saturation energy and 
mode depend on the Rossby (Ro) number. For the hrst set of models either 
oscillatory (with ~ 2 yr period) or steady dynamo solutions are obtained. The 
models in the second set naturally develop a tachocline which, in turn, leads to the 
generation of strong mean magnetic held. Since the held is also deposited into the 
stable deeper layer, its evolutionary time-scale is much longer than in the models 
without a tachocline. Surprisingly, the magnetic held in the upper turbulent 
convection zone evolves in the same time scale as the deep held. These models 
result in either an oscillatory dynamo with ~ 30 yr period or in a steady dynamo 
depending on Ro. In terms of the mean-held dynamo coefficients computed using 
FOSA, the held evolution in the oscillatory models without a tachocline seems to 
be consistent with dynamo waves propagating according to the Parker-Yoshimura 
sign rule. In the models with tachoclines the dynamics is more complex involving 
other transport mechanisms as well as tachocline instabilities. 

Subject headings: Sun: interior — Sun: rotation — Sun: magnetic held — stars: 
interior — stars: activity — stars: magnetic held 


Introduction 


One of the most challenging questions in astrophysics is the origin of the stellar mag¬ 
netism. More specihcally, we do not know yet where inside the stars the large-scale magnetic 
helds are generated and sustained. The best astrophysics laboratory to study stellar mag¬ 
netism is our host star, the Sun. Its magnetic held at photospheric levels can be observed 
in detail, and is depicted in the well known time-latitude “butterhy” diagram. It summa¬ 
rizes the cyclic properties of the surface evolution of the held, its migration patterns, the 
periodicity of magnetic activity and polarity reversals. The morphology and distribution of 
the held in deeper layers is, unfortunately, still evasive to any kind of observation. In other 
solar-like stars (in the main sequence and with spectral types from F to M) the observations 
of magnetic held are less detailed. Nevertheless, cyclic magnetic activ ity, with cycle-period s 
between 5 and 25 years, has been observed in a good sample of stars f Ralinnas et al.11199511 . 
In the same work, the authors reported also observations of hat, non-oscillatory, magnetic 
activity. A general trend is that the strength o f the m agnetic held increases with the rota¬ 
tion rate of the star. More recently, iPetit et al.l (1200811 have been able to infer the topology 
of the surface magnetic held of four stars which could be classihed as solar twins (with 
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stellar parameters close to those of the Sun). Their results indicate that the toroidal held 
strength is proportional to the rotation rate but the poloidal component anti-correlates with 
it. All these observations allow us to look at the stellar magnetism from a broad perspective 
which can add new constraints in the underlying processes and ultimately lead to a better 
understanding of it. 


From the theoretical point of view, the generation and evolution of large-scale magnetic 
helds in stars land other cosmic objects) has been studied using the mean-held dynamo the¬ 
ory (see [Brandenburg fc Subramaniaii 2005 . for a complete review). The mean-held induc¬ 
tion equation implies that the dominant term generating an azimuthal (toroidal) magnetic 
held depends on gradients of the angular velocity, {Bp ■ V)f2, where is the mean angular 
velocity, and Bp is a pre-existing large-scale poloidal held. The generation of the poloidal 
held relies on the turbulent helicities, the so-called a-ehect, a non-dihusive contribution of 
the electromotive force whose dependence on the convective turbulence and rotation is still 
uncertain. Under appropriate conditions, the evolutio n of the magnetic held in this the ory is 
consistent with dynamo waves that, according to the iParkerl (1195511 - lYoshimural (1197511 sign 
rule, propagate in the direction of: 


s = aVfl X Brf 


( 1 ) 


where is the unit vector in the azimuthal direction. The validity of this rule has been 
exten sively tested in two-dimensional kinematic mean-held models (see, e.g. 


Charbonnean 


201011 . For the solar rotation inferred from helioseismology (ISchou et ahl 1199811 . the signs of 
the Vfl components are well known. The sign of a estimated for granular convection is pos¬ 
itive (negative) in the Northern (Southern) hemisphere. These properties favour a dynamo 
process distributed over the entire conv ection zone with th e photospheric evolution sh aped 
by the near-surface shear layer (NSSL) (jBrandenburdl2005l: iPipin &: Kosovichevi 1201111 . Al¬ 
ternative, hux-transport dynamo models in which the surface evolution is shaped by the 
dynamics at the tachocline (and magnetic buoyancy) require coherent meridional motions 
directed equatorward at, or slightly below the tachocline. 

Different observational techniques have reveale d a well-organized, though fl uctuating. 


poleward meridional flow at the solar surface (e.g., IUlrichll2010l: iHathawavi 1201211 with am¬ 
plitudes of about 15 ms“^. However, recent helioseism ology inferences have determined th e 


existence of two or more cells in the radial direction (jZhao et al.ll2013l: ISchad et al.l 1201311 . 
Thus, the existence of a coherent equatorial flow, able to transport the magnetic flux from 
the poles to the latitudes of sunspot activity, is unlikely. 

The dynamo wave propagation sign rule rema ins valid also for mea n-held models con¬ 
sidering the dynamical evolution of the a-effect (iGuerrero et al.l 1201011 . In these models 
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the kinetic component of the a-effect changes in time dne to the backreaction of the mag¬ 
netic field on the plasma motion. The magnetic a-effect depends on the distribntion of 
the small-scale cnrren t helicity, and is a conseqnence of the magnetic helicity conservation 
flPouquet et ahlll976l) . I t is regularly a s snmed that it dimin i shes the indnctiv e action of 
the kinetic helicity (e.g., ISur et ahl 120071: iGuerrero et al.ll2010l: iMitra et ahll201ll) . For high 
magnetic Reynold numbers, this contribution can lead to the total suppression of the dy¬ 
namo action, the so-called catastrophic quenching, in the case when ther e is no effective 
mechanism to remove the sm all-scale current helicity fj Guerrero et ahll2010l) . Recent works 
flVishniac &: Shapovalovll2014ll . and the results that will be presented here, indicate that this 
is not necessarily the only case, and that the magnetic co ntribution to th e a-effect could be 
a source of magnetic field. A similar idea is suggested bv lRonannol (120131) . 


Several kinematic or dynamic mean-held moc 


serve d surface features of the solar magnetic held flDikpati fc Gilman 


els have qualitatively reproduced the ob- 


2001 


2004; Guerrero fc de Gouveia Da l Ping 20081: Pipin fc Kosovichev 201 ll. 


Ghatteriee et ah 


201311 . Based on 


these models, Jouve et ah r fliinoIhlpiDii] fl2014ll have also studied dynamos for solar-like stars 


with the aim of reproducing observed trends for magnetic held strenghts and cycle periods. 
Unfortunately, parameters like helicity and/or turbulent dihusion, in general, are poorly 
determined, and thus the physics of the solar and stellar dynamos cannot be described un¬ 
ambiguously by the mean-held models. An alternative approach is provided by global 3D 
MHD simulations. In this class of modeling, even though the parametric regime is still 
far from the conditions of the stellar interiors, the physics is self-consistently described and 
provides an important insight into the turbulent dynamics and dynamo. 

For instance, systematic hydrodynamic (HD) studies of rotating turbulent convection, 
for the Sun and solar-like stars, have found that the Rossby number, Ro = Mrms/2DoT 


Mrms is the typical turbulent velocity), characterizes 

he resulting mean hows ( 

Gilman 

1976; 

Glatzmaier & Gilman 

19821 Stehen & Frevtae 2007; 

Matt et al. 

2011; 

Kanvla et al. 

2011; 

Guerrero et al. 2013; 

Gastine et al. 2014: Featherstone & Miesch 

2015 

). For large values 


of Ro the models result in an anti-solar diherential rotation (with faster rotating poles 
and a slower rotating equatorial region), and a meridional circulation with single cells per 
hemisphere. In the models with small values of Ro the diherential rotation is solar-like, 
and the meridional how consists of several circulation cells per hemisphere. The transition 
between these regimes is sharp, occurring in a narrow range of Ro. Several properties of the 
observed solar diherential rotation, like the latitudinal dependence of the angular velocity, 
the tachocline 


(IGuerrero et al. 


Ghizaru et ah 


20131: 


boioh. an d a NSSL, although so far only at high latitudes 


Hotta et al 


201511 . have been successfully reproduced. 
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Magnetohydrodyna mic (MHD) globa l dynamo models, which have been developed for 
more than a decade fe.g.. lBrun et al.ll2004l) . have provided a wide spectrum of results, but the 
understanding of the dependence of the mean-flows or magnetic helds on stellar parameters 
is still incomplete. Furthermore, to our knowledge, no MHD model so far has been able to 
reproduce the solar differential rotation. 


MHD simulations of fast-rotating stars by iBrown et al.l (120101) have been able to obtain 
large scale dynamo action. Their results reyealed a steady magnetic field organized in the 
shape of torus, or wreaths, around the equator and with opposite pola rity accros the hemi¬ 
spheres. Non-oscillatory dynamo solutions have also been obtained by ISimitey et ahl (120151) 


Oscill atory dynamo solutions in the global MHD models were first obtained by lGhizaru et ah 


(120101) using the EULAG-MHD code with an implicit sub-grid scale (SGS) formulation. Since 


then other g roups have been able t o simulate magnetic cycles with the use of SGS turbu¬ 


lent models (lAugustson et ahl 120131) . or also in higher-resolution simulations by considering 


only a fraction (a wedge) of the star (IKapyla et al.ll2012l: IWarnecke et ahl 120141) . Although 

perl 


ods obtained in the model of 

Ghizaru et al.l 

Augustson et ah 

2013; 

Warnecke et ah 

2014) 


in this paper, this difference may be explained by t he absence of a ta c hoclin e and a radiatiye 


stable layer in the models of iKapyla et al.l (120121) : lAugustson et al.l (l2013l) : IWarnecke et ah 
(120141) . By consideri ng a forcing function to impose a constant angular yelocity at the bot¬ 


tom of the domain, (iBrowning et ahl 120061) analyzed the effects of a tachocline in convective 
dynamo simulations. Although they found strong toroi dal magnetic held de yelopment at 
this imposed shear layer, they observed no field reversals. iMasada et al.l (120131) explored the 
effects of penetrative convection in spherical dynamo simulations. Due to dissipative effects 
the tachocline obtained in their model is not well defined. However, they obtained stronger 
and cyclic large-scale toroidal magnetic field, demonstrating the importance of the stable 
layer in the storage of the magnetic field. 


In this paper we compare 3D global MHD convective dynamo simulations with different 
Rossby numbers for models with and without the tachocline. Our goal is to compare the 
global properties of the mean-flows and the resultant magnetic actiyity, and contrast them 
with observational signatures. Of particular interest are the inductive and diffusive terms 
resulting from the differential rotation and the collective effects of turbulence. Altogether, 
this will allow us to elucidate on the importance of the tachocline in solar and stellar dynamos. 
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2. The model 


We adopt a full spherical shell, 0 < 0 < 27r, 0 < 6 < ir, the radial domain has its bottom 
boundary at ri, = O.bli?© for the models that develop a tachocline, a nd Vh = 0.72.Rm for th e 
models without it, and the upp er boundary at n = ^ 6 i? 0 . Unlike iGuerrero et al.l (120131) . 
where the anelas tic equations oflLipps &: Hemlerl (119821) were employed, here we solve their 
MHD extension ( Ghizaru et ah 201o[ 


V ■ (p^u) = 0, 


Du 

'm 


2f2 X u = —V 


0 ' 1 
+ 9 -^ H- 

f^Ops 




DO' 


0 / 


Dt 


= -u ■ V 0 , - 


r 


DB 


= {B ■V)u- B{V -u) 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


where D/Dt = d/dt + u ■'V is the total time derivative, u is the velocity held in a rotating 
frame with = f 2 o(cos 6 ', — sin 6 ', 0 ), p' is a pressure perturbation variable that accounts 
for both the gas and magnetic pressure, B is the magnetic held, and 0' i s the potential 
temp erature perturbation with respect to an ambient state 0e (see §3 of iGuerrero et ah 


20131 . for a comprehensive discussion). Furthermore, ps and 0* are the density and potential 
temperature of the reference state which is chosen to be isentropic (i.e., 0 ^ = const) and 
in hydrostatic equilibrium; g = GM/r'^ej. is the gravity acceleration, G and M are the 
gravitational constant and the stellar mass, respectively, and po is the magnetic permeability. 
The potential temperature, 0, is related to the specihc entropy: s = Cpln0 + const. 

The term 0'/r represents the balancing action of the turbulent Reynolds heat hux 
responsible for maintaining t he steady axisy mmetric solution of the stellar structure (see 
section 1.2 and Annexe B in iGossettd l2015l . for details). In this w ork the parame t ers o f 
the ambient state are slightly different from those in equation (7) of Guerrero et ah ( 20131) 
§3.1. For the models including a tachocline, here we use the polytropic indexes = 2 
and rricz = 1.499978, together with the transition width wt = O.OlSi?©. In the models 
without a tachocline, the polytropic index is constant rricz = 1.499985. These values result 
in convective motions with similar Ro for both types of models. The relaxation time of the 
potential temperature perturbation for all models is r = 1.036 x 10® s ( 3.3 yr). 


The equations are solved numericallv using the EULAG-MHD code ( 

Ghizaru et ah 

2010 ; 

Racine et ah 

2011 ; 

Smolarkiewicz & Charbonneau 

2013; 

Guerrero et ah 

2013h. a spin- 

-off of 


the hyc 
search 


rodynamical mod el EULAG predominantly used in atmospheric and climate re- 
Prusa et al.N20081) . The time evolution is calculated using a special semi-implicit 
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method based on a high-resolution, non-oscillatory forward-in-time advectio n scheme MP- 
DAT A (Multidimensional Positive Definite Advection Transport Algorithm; ISmolarkiewicz 
20061) . The truncation terms in MPDATA evince viscosit y comparable to the explicit 


SGS viscosity used in large-eddy simulation 


Domaradzki et al.l 120031: iMareolin et ah 


LES) models ( Elliott fc Smolarkiewiczl 12002 


2006). Thus, the results of MPDATA are often 


interpreted as implicit LES or ILES (ISmolarkiewicz fc Margolinll2007|) . 


For the velocity field we use impermeable, stress-free conditions at the top and bottom 
surfaces of the shell; whereas the magnetic field is assumed to be radial at these boundaries. 
Finally, for the thermal boundary condition we consider zero divergence of the convective 
flux at the bottom and zero flux at the top surface. 


3. Results 

Table 1: Simulation parameters and outputs: Ra* = here Se is the specihc entropy 

of the ambient state, Urms is the volume averaged rms velocity (in m/s) in the unstable layer, 
Ma = Urms/c* is the Mach number, with c* = ■>/ iRT* |r=o. 85 i ?0 being the sound speed at the 
middle of the unstable layer, Ro = and xq = (Deq — r2p)/r2o) where Deq = G(i?Q,0°), 
Dp = D(R 0 , 60°) are the surface rotation rates at 0 and 60 degrees in latitude, and Dq is 
the frame rotation rate expresed in terms of the solar rotation rate Dq. The growth rate of 
the magnetic field. A, is given in Tesla yr“^. The kinetic and magnetic energy densities, in 
Jm“^ are: cr = 'PgU^/2, eu = B /2/io, = B^/2^q and Cp = {B^ + Bg)/2no. Finally, the 
full cycle period, Tm, is expressed in years. Models starting with the letter CZ consider the 
convection zone only, while models starting with RC include both, radiative and convective 
zones. The number of grid mesh points is Nr = 47, Ng = 64 and Ntp = 128 for CZ models, 
and Nr = 64 for RC models. 


Model 

Do 

Ro 

Ra* 

^rms 

Ma(10-4) 

Xo. 

A 



Cp/CK 

7m 

CZOl 

2Do 

0.030 

1.29 

55.12 

3.85 

0.09 

2.26 

0.098 

0.080 

0.018 

2.26 

CZ02 

D© 

0.067 

5.16 

60.91 

4.25 

0.17 

0.65 

0.010 

0.009 

0.001 

2.21 

CZ03 

O 0/2 

0.150 

20.65 

67.18 

4.69 

0.18 

0.56 

2•10-4 

2 • 10-4 

1 •10-5 

- 

RCOl 

2Do 

0.033 

1.36 

60.80 

4.10 

0.07 

0.84 

0.249 

0.042 

0.206 

- 

RC02 

D© 

0.069 

5.45 

62.27 

4.30 

0.05 

0.85 

0.184 

0.163 

0.020 

34.5 

RC03 

D©/2 

0.161 

21.80 

72.88 

5.04 

0.28 

0.06 

0.004 

0.004 

3•10-4 

- 


We have performed two sets of simulations: 1) only for the unstable stratihed convection 




































zone (CZ models), and 2) for the convection zone with a convectively stable radiative zone 
at the bottom of the domain (RC models). The hrst set does not support the formation of a 
strong radial shear at the base of the convection zone (the tachocline), and thus excludes it as 
a source of magnetic held. The second set naturally leads to the development of a tachocline 
and therefore to the {B ■ V)f2 source of the toroidal held acting in the tachocline. The input 
and resulting parameters of the models for three diherent rotation rates, flo (corresponding 
to diherent Rossby numbers )0 are summarized in Table [H 


^Since we change the value of the Ro only by changing the value of the frame rotation rate, use both 
terms interchangeably to express the dependence of the results with flo- 



9 


3.1. Dynamo models without tachocline 


In Figure [T] we present, from left to right, the meridional prohles of the differential 
rotation, meridional circulation, and snapshots of the vertical velocity, Ur-, and the toroidal 
magnetic held resulting from the models a) CZOl, b) CZ02 and c) CZ03. The differential 
rotation and meridional circulation correspond to averages over longitude and time (~ 3 
years). The instantaneous orthographic projections of Ur and allow us to distinguish the 
character of the convective how and the distribution of the toroidal held for each model. 
Unlike our previous purely hydrodynamic simulations, in which we found a c ritical Ro value 


divid ing the rotation prohles into solar-like and anti-solar rotation types fj Guerrero et ah 


20131 ). in the present simulations with magnetic held all the models exhibit the solar-like 
diherential rotation, irrespective of the Rossby number. This proper ty is related to the 
inhuence of the dynamo-generated magnetic held on the huid dynamics. iKarak et al.l (120151) 
have found that for the MHD case the transition occurs at higher values of Ro. For instance, 
model CZOl shows iso-rotation contours mainly aligned along the rotation axis (see Fig. [T^). 
In models CZ02 (Fig. [T]d) and CZ03 (Fig. [T}:;) the rotation prohle exhibits clear conical-shape 
contours that resem ble the rotation of the solar convection zone inferred by helioseismology 
(jSchou et al.lll998f) . During early stages of the models evolution, when the magnetic held 
inhuence is negligible, the rotation contours are aligned along the rotation axis (the so- 
called Taylor-Proudman balance), and then transformed into the conical-shape contours. 
This indicates that both Reynolds and Maxwell stresses contribute to the distribution of 
angular momentum. Later, it is remarkable that even model CZ03, in which the hnal 
magnetic energy is 10^ times smaller than the kinetic energy, the rotation law departs from 
the Taylor-Proudman balance. A detailed discussion of the angular momentum balance is 
out of the scope of this paper and will be addressed in a subsequent paper. 
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a) CZOl 


b) CZ02 




c) CZ03 




Fig. 1.— From left to right: mean profiles of the rotation rate and meridional circulation, 
and snapshots at r = 0.95Rq of the vertical velocity, Ur, and the toroidal held, for 
the simulated models a) CZOl, b) CZ02, and c) CZ03. In the meridional circulation 
panels continuous (dashed) lines correspond to clockwise (counterclockwise) circulation. The 
background hlled contours show the latitudinal velocity, ug. The prohles of differential 
rotation and meridional circulation correspond to mean azimuthal values averaged over ~ 3 
years during the steady state phase of the simulation. 
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Another remarkable feature of models CZ02 and CZ03 is the natural development of 
a well-defined NSSL. The radial shear in this layer is negative and extends fr om the equator 


to the poles as observed in the Sun. In agreement with the results found in iGuerrero et ah 


(120131) ■ the NSSL arises from the appropriate choice of the ambient state, ©e, and the 
relaxation time of the potential temp erature perturbations , r. As mentioned above, these 
quantities differ from the ones used in I Guerrero et ahl (120131) in such a way that Ro becomes 
larger at the top of the domain (i.e., in a thin layer close to the surface where the buoyancy 
force becomes much stronger than the Goriolis force). Although the number of mesh points 
is not sufficient to resolve super-granulation scales which may be important for the formation 
of the solar NSSL, our model succeeded in reproducing similar effects in the resolved scales. 
Furthermore, the latitudinal angular velocity gradient in both models corresponds well to 
the solar rotation, i.e., xn = (f^eq — f^p)/f^o = 0.18, where Geq is the equatorial angular 
velocity, Gp is the angular velocity at 60° latitude (both quantities are computed from the 
temporally and azimuthally averaged profile of G), and Gq is the angular velocity of the 
reference frame. 


As far as the meridional circulation is concerned, models GZOl and GZ02 show a 
multicellular pattern (second column of Fig. [T^ and b), while model CZ03 shows a dominant 
(counter) clockwise cell in the (Northern) Southern hemisphere with a latitudinal velocity, 
Mg ~ 12 m s“^ (about twice the value of ue in models GZOl and GZ02). This cell is located 
near the base of the convection zone. Another cell, with a smaller velocity and opposite 
circulation direction appears in the subsurface layer (second column of Fig. [Tt). 


An oscillatory large-scale dynamo action is observed in models GZOl and CZ02 as 
can be seen in Fig. [2^,b. In both cases the toroidal magnetic field is generally symmetric 
across the equator and reverses with a cycle period of ~ 2 years (Tabled]). For the model 
GZ02 we found a secondary ~ 6-years cycle modulating the amplitude of the magnetic 
field. These properties are at odd with the solar 22-year cycle but in agreement wi th other 


global models without tachocline (e.g.. lAugustson et al.ll2015l: IWarnecke et al.ll2014]) . In the 


convection zone the evolution of the magnetic fields is consistent with a pattern of dynamo 
waves propagating from the bottom to the top of the domain (Fig. [2^,b). At the surface 
the field forms a ±50° belt around the equator. The latitudinal migration of the field is 
slightly poleward. The long term evolution of the simulations GZOl and GZ02 indicates 
extended periods of minimal or maximal activity (see for instance the extended minimum 
in the Northern hemisphere of the model GZOl, Fig. [2^, between 20 and 30 years). The 
amplitude of the dynamo-generated magnetic field depends on the rotation rate, reaching 
~ 10% and ~ 1% of the kinetic energy of the system, for the models GZOl and GZ02, 
respectively (see Table [1]). In model GZ03 (Fig. the magnetic field grows exponentially 
but saturates at ~ 0.02% of the kinetic energy, (see Tabled]). The magnetic field in this case 
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Fig. 2.— Time-latitude (left) and radius-latitude (right) diagrams of the toroidal magnetic 
held, Bfj,, for models a) CZOl, b) CZ02 and c) CZ03. On the left, the contours show B^, in 
Tesla, at r = O.95i?0. On the right, the contours are taken at 30° latitude. Only a fraction 
of time of the simulated statistically steady state is shown in each plot. 
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is nearly steady. Its time evolution shows periods when the held amplitude abruptly decays, 
stays in a minimum state for some years, and then quickly grows again to its saturation 
amplitude without polarity reversals (Fig. [2t). As indicated in Table [H the growth rate of 
the dynamo is proportional to the rotation rate. 
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c) a.Q/?oV^to. cz 

2000 

1000 

0 

-1000 

0.61 0.70 0.78 0.87 0.96 0.61 0.70 0.78 0.87 0.96 

7-/7?0 r/i?0 



d) a,ni?oV77„, Kc 



e) r 'dgQRj'/rj^o, CZ 



90-e 


f) r-'a,QfteVdto. RC 



90-6 


Fig. 3.— Magnetic field source terms (see Eqs. 8-9) for models CZ (left panels), and RC (right). 
Panels a) and b) show the FOSA estimation for aRoi/rjto- Dotted and dot-dashed lines represent 
the kinetic and magnetic contribution, respectively. The continuous line is the resultant, non- 
dimensional, total a-effect. The profiles of a are latitudinal averages in the northern hemisphere. 
Middle and bottom rows show drRlR%/r]iQ and r~^deQ.R^/r]iQ, respectively. In panels c) and d) 
continuous (dashed) lines indicate latitudinal averages at lower, 20° — 40°, (higher, 60° — 80°), 
latitudes. (Note that in panel d) the value of dr^R^ has been divided by 2.) In panels e) and f) 
continuous (dashed) lines are vertical averages in the top, 0.9Rq < r < 0.95 Rq (bottom, O.72i?0 < 
r < O. 77 R 0 ), of the domain, mn is a radial average of defined in Eq. (19]), for r > 0.72720. 
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Table 2: Dynamo coefficients compnted from the simnlations resnlts. The dynamo numbers 
are computed with Eqs. ([8]) and ([9]). The 0 index in Ca and Cn refers to the maximum 
value of each quantity, in rjt, to its average over the CZ (in m^s“^). The dynamo numbers 
Dr = CcioCqq and Dq = CaoC^Q, provide an indication of the dynamo efficiency when 
considering the radial and the latitudinal shear, respectively. Finally max(ile) 

and max{Br) are the maximum absolute values of the large-scale components of the field at 
r = O.95i?0 (the values in brackets are measured at r = O.72i?0) in units of Tesla. 


Model 

mo(10®) 

/^k 

^aO 

/^m 

^aO 

CaQ 

D^o(lO^) 


Dr 

De 

niax(i3,^) 

max(i3g) 

max(Br) 

CZOl 

0.99 

5.01 

1.84 

4.47 

0.37 

27.28 

1.65 

122.90 

0.60 

0.18 

0.03 

CZ02 

1.09 

5.38 

1.28 

5.00 

1.62 

14.75 

8.10 

73.75 

0.61 

0.13 

0.06 

CZ03 

1.28 

2.72 

0.08 

2,72 

1.69 

9.71 

4.59 

26.41 

0.15 

0.05 

0.02 

RCOl 

1.08 

3.03 

2.06 

2.75 

0.81 

19.65 

2.22 

54.03 

0.52(0.19) 

0.10(0.80) 

0.04(0.26) 

RC02 

1.12 

2.80 

1.52 

2.72 

1.67 

13.33 

4.54 

36.25 

0.22(0.94) 

0.07(0.28) 

0.03(0.24) 

RC03 

1.25 

3.56 

0.64 

3.52 

4.98 

21.51 

17.53 

75.71 

0.20(0.71) 

0.04(0.37) 

0.03(0.20) 


The spatio-temporal evolution of the magnetic held in the models can be interpreted in 
terms of its main sources according to the mean-held dynamo theory, i.e., the a and D ehects. 
The prohl e of a is estimated by using the hrs t-order smoothing approximation (FOSA, see 
for details iBrandenburg fc SubramanianI 120051 ) as follows: 




Du 


+ -hc 
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( 6 ) 


where Tc = Hp/urms is the turnover time of the convection, = d\np/dr is the density 
length scale, and 


hk = u' ■'V X u' (7) 
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hoPs 


are the small-scale (note the prime over u and b) kinetic and current helicities, respectively 
(the overbars denote azimuthal and temporal averages over 3 years). As it is common in the 
mean-held dynamo theory, we present the source terms them as ratios between the induction 
and dihusion times: 
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where r/to is the FOSA estimation of the turbulent magnetic diffusivity, 

_ 1 2 
VtO 2 A'i^rms 5 


( 9 ) 


averaged over the entire convection zone, O.72i?0 < r < O.96i?0 (the values of rjtQ are 
presented in Table [2]) . 


In the top, middle and bottom rows of Fig. [3] we present the normalized prohles of 
a, drO. and respectively, computed from the simulation results. Fig. |3^ shows the 

radial prohles of the kinetic and magnetic contributions to a presented by dotted and dot- 
dashed lines, respectively, while the resultant a is shown with continuous lines. The vertical 
prohles of the radial shear. Fig. [3t, are shown at two diherent ranges of latitude: a higher 
latitude prohle (dashed lines) was computed averaging over 60° and 80° latitude, and the 
lower latitude prohle (continuous line) was averaged between 20° and 40° latitude. Finally, 
the latitudinal shear was also computed for two diherent ranges of depths. Fig. |3^ . The 
bottom prohle (dashed line) corresponds to an average between O.72i?0 < r < O.77i?0, while 
the top prohle (continuous line) is an average between O.9Oi?0 < r < 0.95 Rq. Presenting 
these quantities as radial prohles, instead of contours in the meridional plane, allows to 
straightforwardly compare the diherent dynamo models. The red, blue and green lines 
correspond to models 01, 02 and 03, respectively. Furthermore, Figure |3] allows us to 
directly compare the CZ and RC models (to be discussed later), shown in the left and right 
panels, respectively. 


Before describing the evolution of the magnetic helds, it is instructive to discuss how 
the source terms change as a function of the Rossby number. Table [2] presents the maximal 
absolute values of the quantities depicted in Fig. [3] (denoted with a 0 index), as well as 
the maximal amplitude of the mean magnetic held components at the surface level. (Some 
values at the boundaries can be very large without really contributing to the dynamo, so we 
have disregarded them.) In general, there are no major diherences between the prohles of 
the kinetic a-term, In all the cases it is positive in the bulk of the convection zone and 
changes sign near the boundaries. The radius at which this change of sign occurs is shallower 
(deeper) for the faster (slower) rotating model. There is not a clear relation between C^q 
and Ro, it is only noticeable, however, that the faster models, CZOl and CZ02, have 
larger kinetic helicity than the slower model, CZ03. As for the the magnetic a-ehect, 
we should hrst notice that its amplitude is comparable to the kinetic a. However, C^q is 
inversely proportional to Ro. This is expected since the models with faster rotation develop 
stronger magnetic fields. Correspondingly, the total a-effect is inversely proportional to the 
Rossby number. For the models CZOl and CZ02 the signs of and C™ in the middle of 
the convection zone are opposite. At the top of the domain the two terms have the same 
sign, and both contribute positively to the generation of magnetic held. In model CZ03 
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the contribution of C™ is negligible, thus, only the kinetic helicity contributes to the held 
generation. 

Regarding the rotational shear terms, we notice that for larger values of Ro, the surface 
radial shear is stronger and the latitudinal shear is weaker. Note that the latitudinal prohle 
of Cq (Fig. |3^) for model CZOl shows a strong shear at the equator and the poles. This shear 
is reduced in models CZ02 and CZ03 as the vertical motions become more important. These 
fast vertical hows, in turn, are less ahected by the Coriolis force, resulting in a strong radial 
shear in the NSSL. The surface values of the magnetic held are also inversely proportional 
to Ro, possibly indicating that their main source is the latitudinal shear. 

The evolution of models CZOl and CZ02, presented in the butterhy diagrams of Fig. 
|2^,b, is similar. They exhibit branches of magnetic held migrating in latitude from the 
equator towards the poles. This pattern agrees with the Parker-Yoshimura sign rule, Eq. 
(11]), for positive values of adr^ in almost the entire convection zone (see continuous red 
and blue lines in Fig. |3^, c). In the vertical direction the propagation of the magnetic 
held in the middle of the CZ also agrees with Eq. ([T]) for adgQ > 0. Because of the large 
values of Cq in model CZOl, the branches of toroidal held show the large magnetic held 
strength from the bottom of the convection zone to the top. In model CZ02, with a smooth 
latitudinal shear the stronger vertical branches of start at a depth ~ O.SSRq (compare 
the right panels of Fig. |2^ and b). At the upper radial levels of model CZ02, the term 
adoQ becomes negative, and this should change the direction of migration. However, at 
the latitude depicted in Fig. [2l only a slight change in the tilt of the branches is evident. 
The morphology of the azimuthally averaged magnetic held of model CZ02, is depicted in 
Fig. H^-d. The hlled contours show the positive (yellow) and negative (blue) toroidal held 
strength. The continuous (dashed) lines represent the clockwise (counterclockwise) poloidal 
held lines. The four snapshots cover one magnetic held reversal. In the same plot it can 
be observed that the poloidal held follows a wave-like evolution pattern similar to B^. The 
pattern starts to develop at the bottom of the domain at lower latitudes. As the cycle evolves 
it grows and expands over the convection zone with the more intense regions (the innermost 
contours) migrating poleward and upward. In this case, however, the held is antisymmetric 
with respect to the equator. 
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a) t= 24.7 yr b) t= 25.1 yr 



c) t— 25.5 yr d) t- 25.9 yr 





Fig. 4.— Snapshots of the mean magnetic helds, (color images), and poloidal magnetic 
held (contour lines) for the model CZ02. Continuous (dashed) lines represent clockwise 
(counterclockwise) magnetic held direction. The time series covers one polarity reversal 
(half dynamo cycle). 
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3.2. Dynamo models with tachocline 


In this section we present three models (RCOl, RC02 and RC03) with the tachocline. 
All of them have the same stratification and differ only by the rotation rate, ffo- For all 
these cases, the magnetic field evolves differently than in the CZOl, CZ02 and CZ03 models. 

It is worth mentioning first that due to the presence of the tachocline, most of the 
magnetic field develops at the base of the convection zone. Thus, the Maxwell stresses 
play an important role in the downwards transport of the angular momentum. In Figure 
[5l panels a, b and c show the differential rotation (left column) and meridional circulation 
(right column) profiles of models RCOl, RC02 and RC03, respectively. Model RCOl has 
the lowest value of Ro, thus it has the strongest influence of the Coriolis force. This force 
tends to homogenize the rotation of the convection zone with that of the radiative zone. 
This happens in a region between ~ 5° and ~ 80° latitude. However, zones with some radial 
and latitudinal shears are developed at the equatorial and polar latitudes. In model RC02, 
both radial and latitudinal differential rotations develop in larger zones. The intermediate 
latitudes are in iso-rotation with the radiative zone. 


The transport of angular momentum to the radiative zone makes the stable layer to 
rotate, on average (cf. Fig. Od), faster than the frame. Hence, the equatorial acceleration 
of the convection zone in the model RC02 is not as pronounced as in the model CZ02. In 
the former case the angular velocity is around 440nHz, while in the latter case it reaches 
500nHz. Consequently, the latitude of iso-rotation with the frame is moved up to ~ 70°. 
In the model CZ02 it is located at ~ 50° (the yellow filled contours in Figs. [T]d and [5 ]d 
correspond approximately to the frame rotation rate, hlo)- 


Besides affecting the average profile of the differential rotation (with respect to purely 
hydrodynamic models) the magnetic feedback on the flow generates torsional oscillations in 
the models where the dynamo is periodic (CZ02 and RC02). The latitudinal morphology of 
the simulated oscillations resemble the solar observations, but the amplitude of the observed 
oscillations is a fe w times sma l ler. T his might be a consequence of a rather large a-effect 
such as found bv ICovas et al.l (120041) . The origin of the torsional oscillations in our sim¬ 


ulations seems to be due to a modulation of the latitudinal angular momentum transport 
mediated by the meridional circulation and the magnetic torque (at equatorial latitudes). 
A comprehensive discussion of these results is beyo nd the scope of the current work and is 
presented in a separate paper flCuerrero et al.l 120 161) . 


A NSSL, more pronounced at higher latitudes, is also observed in this model (RCOl). In 
both RCOl and RC02 cases the contours of rotation form conical shapes in the convection 
zone (see the vertical dotted lines in Fig. [5] as an eye guide). Finally, in model RC03 there 
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is a well defined latitudinal differential rotation in the bulk of the convection zone, and also 
a NSSL is formed. The region of iso-rotation of the convection zone with the radiative core 
occurs at higher latitudes (~ 70° latitude). For this reason the radial shear at the equator, in 
the tacocline, is stronger in model RC03 than in models RCOl and RC02. The contours of 
model RC03 appear cylindrical which is surprising because in this model the Coriolis force 
influence is smaller than in the other two cases. The rightmost columns of Fig. [5] show the 
convective structure, represented by the vertical velocity and the distribution of the toroidal 
field at r = O.QSRq. 
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a) RCOl 



b) RC02 



c) RC03 




Fig. 5.— Same as Fig. [T]for the models a) RCOl, b) RC02, and c) RC03. For these cases 
the temporal average is over ~ 10 years during the steady state phase of the simulation. 
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Fig. 6.— Same than Fig. [2] for the models a) RCOl, b) RC02 and c) RC03. The color 
scales show B^'m. Tesla. 


The meridional circnlation is mnlticellnlar in all onr models with the tachocline. How¬ 
ever, as expected, model RC03 has larger meridional velocities than the other models and ex¬ 
hibits a dominant connterclockwise (clockwise) cell in the Northern (Sonthern) hemispheres. 
Noteworthy, in all these models a low amplitnde poleward flow develops in the npper part of 
the convection zone at latitndes > 30°. This can be noticed in the second colnmn of Fig. [5] 
in the color hlled contours of the latitudinal velocity, uq^ as the blue and yellow regions in the 
Northern and Southern hemispheres, respectively. The formation of this flow is likely due 
to the gyroscopic pumping mechanism because of the negative gradient of angular velocity 
f Miesch fc Hindman 2011 1. 


The third column of Fig. [5] shows that the smaller convective structures in our simu¬ 
lations are formed at higher latitudes. In the equatorial band the structures are elongated, 
resembling the so-called banana cells. Note that the vertical domain of our simulations 
reaches up to 0.96Rq, and small-scale structures are not resolved. We surmise that this 
is the reason why NSSL and latitudinal meridional circulation only appear above ~ 30° 
latitude. Nevertheless, the meridional flow structure qualitatively agrees with recent helio¬ 


seismology invers ions made using the Solar Dynamics Observatory data (jZhao et al.l 12013 
Schad et ah 2013). 


The magnetic held in these models evolves on longer timescales than in the CZ models 
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(see ^3.3p . The fast rotating model, RCOl, develops a magnetic field that oscillates in 
amplitude but does not show clear polarity reversals. The oscillation period is ~ 10 yr (not 
included in Tabled] because it does not correspond to magnetic held reversals). The topology 
of the held consists of wreaths of toroidal held of opposite polarity across the equator. 

Model RC02 presents magnetic cycles with a full period of ~ 30yr. Unlike the solar 
magnetic held, in this model the toroidal (poloidal) component of the magnetic held is 
symmetric (anti-symmetric) relative to the equator (see Figs. O) and [7|). Most of the 
magnetic energy in this case is in the toroidal component of the held. In the simulation with 
the slowest rotation rate (model RC03), the magnetic held is nearly steady (as in CZ03) 
with a larger concentration of in a narrow region at the base of the convection zone at 
the equatorial latitudes. The held at the surface has the same polarity as at the bottom and 
is concentrated at higher latitudes. 

An important disparity between the models RC and their non-tachocline counterparts, 
models CZ, is the ratio between their volume averaged magnetic and kinetic energies as 
can be seen in Table [H This diherence is absent in the surface layers where both sets of 
models show magnetic helds of similar magnitude (in fact, models CZ have slightly higher 
values, see Table |2|). In both, CZ and RC cases, the amplitude of the surface held correlates 
directly with the rotation rate. Nevertheless, because of the radial shear at the tachocline, 
in the RC models the maximum mean-held in that region reaches up to 1 T (10'^ Gauss), 
as evident in the right panels of Fig. Eland Table |2] (values in parenthesis). 

As for the turbulent coefficients in these models, there is a clear correlation between the 
turbulent dihusion coefficient, ritQ, and the rotation rate, indicating that rotation quenches 
turbulent dihusion. However, while Ro changes by a factor of ~ 5 between the faster and 
slower models, r]to changes only by a factor of ~ 1.15. Similar to the models CZ, there is no 
correlation between the kinetic a ehect and the rotation. Figure Eb shows that the radial 
prohles of are roughly the same for the three models. They have small positive values 
at the radiative zone, reach a positive maximum in the middle of the convection zone and 
change sign near the upper boundary. On the other hand, changes proportionally to the 
rotation rate. The sign of the C™ at the tachocline is positive, and it is the dominant term 
in the total a-effect in that region. In the lower/middle part of the convection zone the sign 
of C™ is opposite to C^, but its amplitude is small. For models RCOl and RC02 in the 
upper convection zone C™ is negative and, because of the low density, reaches maximum 
values. Similarly to models CZ, both and C™ contribute to the held generation. For the 
slowest rotating model, RC03, the contribution of is unimportant. 

Models RC develop strong radial shear at the interface between the radiative and con¬ 
vective layers (Fig. Eli)- This radial shear anti-correlates with the rotation rate, i.e., the 
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slower the rotation rate, the largest the shear at the tachocline. Similar to the solar case, the 
radial shear at the tachocline is negative at the poles and positive at the equator. According 
to helioseismology the Sun exhibits the strongest shear at the poles, whereas all the models in 
this section have the strongest shear at the equator. The prohles of drfl are fairly flat in the 
convection zone and negative in the NSSL at high latitudes (Fig. |3ll). As for the latitudinal 
differential rotation, the fast rotating model shows a peak at low latitudes (between 15° and 
40° latitude) followed for a flat, zero shear, region in middle latitudes and larger gradients 
of in latitudes above 60° (see red line in Fig. [3]F). Models RC02 and RC03 have prohles 
which are comparable at lower and middle latitudes. However, near the poles, the prohle of 
model RC03 shows a stronger latitudinal shear. 

The morphology of the steady magnetic held in models RCOl and RC03 can be described 
in terms of the distribution of its source terms. For model RCOl, at the tachocline, both 
sources of magnetic held are small, the dynamo results in a steady held concentrated at the 
equator where peaks (red dashed line in Fig. [3]F). In this region the a-ehect is possibly 
dominant over the H-ehect, yet the poloidal component of the held is ~ 4 times larger 
than the toroidal one. In the upper convection zone, has larger values at the equator, 
which explain why the held is concentrated in a band at ±30° latitude (Fig. [6^). The 
held variations do not lead to polarity reversals, perhaps because the held at the tachocline 
remains steady. For model RC03 the held reaches its saturation after t ~ 340 yr with 
steady magnetic helds. Since in this case d^fl reaches larger values at the poles, the toroidal 
magnetic held is localized at polar latitudes. 
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03 


c) t= 60.5 yr 


d) t= 65.8 yr 



a) t= 50.0 yr 


0,00 




Fig. 7.— Same than Fig. IHfor model RC02. 
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Unlike the oscillatory models without the tachocline, CZOl and CZ02, the magnetic 
held migration in model with the tachocline RC02 cannot be solely explained in terms of 
the Parker-Yoshimura rule. Because of the long cycle-period of activity, other transport 
processes like the meridional circulation or turbulent pumping can affect the magnetic held 
evolution. At the tachocline the migration of the held can be explained by the Parker- 
Yoshimura rule (Eq. [1]). The held developed at the poles migrates equatorward, and the 
held developed at the equator migrates poleward. This migrating process can be seen in 
Fig. [7^-d which, similarly to Fig. HJ shows a half of the dynamo cycle with one magnetic 
held reversal illustrates in four snapshots. These dynamo waves result from a positive a and 
negative (positive) values of dr^l at higher (lower) latitudes. The upward radial migration 
from r ~ O.TSi?© to r ~ O.QOi?© is also in agreement with Eq. ([1]) for adgQ > 0 in this 
region. In the upper convective layer, above r ~ O.9Oi?0, there is a region of strong magnetic 
held, possibly resulting from the large negative radial shear. The latitudinal migration in this 
region, depicted in Fig. [6]d, shows two branches of slighty equatorward migration: one branch 
is located between 0° and ±30° latitude, the other one is above 50° latitude. According to 
Eq. (II]) these branches would require negative values of s; however this does not agree with 
the prohles of a and drfi shown in Fig. 12b and|Hb (blue lines). The most equatorial branch 
of coincides with a clockwise meridional circulation cell with an equatorward how at the 
surface (Fig. [5b). However, it is impossible to assure that this is the case giving the difficulty 
in disentangling the local and non-local ehects contributing to the held migration. 


3.3. What process governs the cycle-period? 


Perhaps the most interesting diherence between the oscillatory dynamo models CZ02 
and RC02, is the timescale of the magnetic cycles (i.e., the period of the magnetic polarity re¬ 
versals). As follows from the mean-held dynamo theory, the magnetic held grows and decays 
according to the values of the indu ctive and the dihusive terms in the induction equation 
ferandenburg fc Subramanian 2005 1. A long standing problem for mean-held dynamo mod¬ 
elers has been to conciliate the theoretically expected value of the turbulent dihusivity, rjt, 
with the observed cycle period. In the solar case a mixing-length theory (MET) estimation of 
rjt gives values ^ 10^ m^s~^ which, in turn, results in dynamo cycles with periods of 2 — 8 yr 


(IGuerrero et al.ll2009fl instead of the observed 22 yr. Non-linearities in the dependence of the 
turbulent dihusivity on the large-scale magnetic held, li ke the so-called p-quenching mecha 


nism , have been explored but did n ot solve the problem flRiidiger et al.lll994 iGuerrero et ah 


20091: iMuhoz-Jaramillo et al.ll201lh . The cycle-period diherence between the two types of 


models found in our simulations provides an opportunity to explore this issue. 
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The radial profile of the turbulent diffusion coefficient, rjt, computed using Eq. (|9]), for 
both models with the solar rotation rate, CZ02 and RC02, is depicted in Fig. |8l From this 
hgure and the values of r^to in Table [U it can be seen that the prohles and values of rjt in the 
convection zone match closely in both models. Both cases correspond to a strongly diffusive 
regime with r]t ~ 10® m^s“^ (in agreement with the MLT estimate). However, model CZ02 
has a ~ 2-yr cycle period, while model RC02 has a full magnetic cycle-period of ~ 30 years. 
Furthermore, in the same hgure, the dashed lines correspond to the turbulent diffusivity 
of the hydrodynamic versions of the CZ02 and RC02 models. We notice that there is no 
important diffusivity quenching due to the presence of the saturated magnetic held except 
in a small fraction of the radius nearby the tachocline where the dihusivity of the magnetic 
model is smaller (see the region between the vertical dotted lines in Fig. [8]). 



r //?0 


Fig. 8.— Vertical prohles of the turbulent magnetic dihusivity, rjt, for model CZOl (blue) 
and RC02 (red). The dashed lines show the prohles of r]t for the hydrodynamic versions of 
models CZ02 and RC02. 


This suggests that the cycle period in model RC02 may be determined by the value of 
rjt in the region where most of the magnetic held is produced and stored (i.e., at and below 
the tachocline), which is mainly convectively stable. (Note, that T]t does not tend to zero in 
the radiative layer but to values around ~ 10® m^s“^. This reveals that the magnetic held 
is inducing turbulent motions in some fraction of the radiative interior.) To verify if this is 
indeed the case, we have computed the ratio between the dihusion times for models RC02 
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and CZ02, as follows: 


i-RC 

i-CZ 


L%c{Vt)cz 

Lcz{Vt)RC 


~ 19 


( 10 ) 


where, Lcz = (0.96 — O.72)i?0 is the length scale of the convective layer of model CZ02, and 
Lrc = (0.77 — O. 68 )i ?0 is the length of the region of model RC02 where most of magnetic 
held is stored (see vertical dotted lines in Figs [3^ and b); {rit)cz and {rit)Rc are the radially 
averaged values of the turbulent diffusivity over the same lengths. Interestingly, this rough 
estimation agrees fairly well with the ratio between the cycle periods of these two types of 
models. 

It is striking that in model RC02 the magnetic held in the turbulent convection zone, 
despite the larger local values of T]t, evolves on the same timescale as the magnetic held 
in the most stable layer with lower dihusivity. If such a non-Iocal mechanism is operating 
in the Sun, it could explain the 22-year period of the solar magnetic cycle. This does not 
necessarily mean, however, that the migration of the sunspots observed at the surface is 
shaped by the meridional circulation at the base of the convection zone as assumed by 
most of the hux-transport dynamo models. On the contrary, the dynamo in our simulations 
operates distributed over the entire convection zone, and several processes could inhuence 
the magnetic held generation and mig r a.tion (e.g., the near- s urface shear rotational shear 
layer, as argued by BrandenbnrgI (2005); Pipin fc Kosovichev ( 2011 1). 


3.4. Tachocline instabilities 


From Fig. [ 8 ] it is evident that turbulent motions must be present in the stably stratihed 
layer. Such motions increase the turbulent dihusivity by more than one order of magnitude 
compared to the pure hydrodynamic model RC02 (red dashed line). The most probable origin 
for this hydromagnetic turbulence in the stable layer is the development of MHD instabilities 
at and below the interface between the convective and the radiative layers. It has been 
suggested that these turbulent motions could be signihcantly helical suc h that, similarly to 
the h ow in the convection zone, they can result in a kinetic a-ehect fe.g.. iDikpati fc Gilman 
200ll) . The analysis performed in 1 13.21 indicates that the kinetic part of the a-ehect in the 
stable layer is small when compared with its values at the convection zone. Nevertheless, 
the magnetic contribution to a in the radiative zone is important, as evidenced in Fig. [3 ]d. 
The development of small-scale helical magnetic structures due to instabilities could explain 
the existence o f this magnetic g-ehect. This, in tu rn, leads to the generation of large-scale 
magnetic held ({Brandenburg fc Subramanian 2005 1. 


This section does not intend to explore in detail the development of all possible tachocline 
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instabili ties in the simulation s presented here. Such analysis has been recently performed al¬ 
ready bv iLawson et al.l (120151) for one of the EULAG-MHD simulations fjPassos fc Charbonneau 
2014J) . In this section we briefly review the nature of such instabilities, try to identify which 
of them could be operating in our simulations based in a simple energy analysis, and discuss 
their contribution to the dynamo mechanism. 

Tachocline instabilities have been extensively investigated over the past decades. Ac 


cordi ng to the energy source they can be divided into four distinct kinds which include flArlt 


20091 1: a) shear-driven instabilities; b) baroclinic instabilities; c) buoyanc y driven instabili 


ties; and d) current-driven instabilities. The shear-driven instabilities fe.g.. lDziembowski fc Kosovichev 


19871: ICharbonneau et al.lll999l: iDikpati &: Gilman! 1200111 are purely hydrodynamic. Due to 
the latitudinal shear at the tachocline, small perturbations can destabilize the flow given 
sufficiently strong shear, e.g., more than ~ 30% difference betwee n the equator and the pole 
for a 2D analysis. The 3D stability analysis of lArlt et al.l (120051 1 has sh own tha t the s olar 
tachocline should be nearly stable, in agreement with previous results bv iGaraudI (1200111 . In 
the simulations presented in 1 13.21 the difference in the angular velocity between equator and 
pole is less than 10% for the models RGOl and RG02, and about 30% for the model RG03. 
It seems unlikely that this kind of instability is developing in these simulations. Baroclinic 
instabilities (b) are more complex in nature. If the transition region from subadiabatic to 
super-adiabatic is in thermal wind balance, horizontal (azimuthal and latit udinal) pertur- 
bation flows can gain some thermal energy and become unstable. Recently iGilmaru (120151) 
has demonstrated that the magnetic held has a stabilizing effect on these growing modes. 
Buoyancy driven instabilities (c) are presumably responsible for the emergence of magnetic 
hux tubes. The buoyancy force of low density parcels of magnetized gas is responsible for the 
emergence. For toroidal held ropes stored in a stable layer, it is generally tru e that the mag¬ 


netic held should surpass the equipartition held strength (see iHughesI 120071. and references 


therein). This is not the case in the simulations above. Strong hux concentration could be 
achieved in the case of very thin hux tubes which would require higher grid resolutions. The 
current-driven instabilities (d) appear whenever there is a large-scale current and their en¬ 


ergy source is the large sca l e magnetic held 


his k ind of instabilities could develop withou t 


rotation fe.g.. lTavlerlll973l: iBonanno fc UrDiijl2012l) or with r otation (jPitts fc Tavlerlll985h . 


or also when the rotation is dihere ntial (e.g.,lArlt et al.ll2005[) . These latter processes, often 
called magneto-shear instabilities ( Gallv et ah 20031) . destabilize diherent conhgurations of 
toroidal magnetic bands either by opening the magnetic held lines (the so-called clamshell 
instability) or by tipping the axis of the magnetic held band. In both cases it is expected the 
development of non-axi-symmetric magnetic held out of the axi-symmetric held lines. This 
instability grows for diherent conhgurations of toroidal magnetic helds and diminishes the 
shear. When, a latitudinal shear prohle and a toroidal held conhguration are left to evolve 
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freely, after the saturation the instability decays fjCallv et ahll2003l) . However, if the system 
is forced, as the solar tachocli ne should be, both the shear and the magnetic held adjust to 
a new equilibrium state 


fe.g. 


Miesch et ahll2007[) . 


All the ingredients for the onset of the magneto-shear instability are present in the 
simulations with radiative zone presented above (RC models). The differential rotation is 
sustained by the Reynolds stresses and, depending on the Rossby number, it induces different 
conhgurations of toroidal helds (which in turn react back re-adjusting the shear). At some 
evolutionary stage of the simulation the conditions for the development of the instability are 
present creating small scale hows and magnetic held in the stable layer. As a consequence, 
a new adjustment should happen until the hnal steady state of the simulation is achieved. 


Following iMiesch et al.l (120071) . we study the time evolution of the diherent energy 
sources of the stable layer to verify the existence, and relevance, of s uch small scale, non - 
axisymmetric, magnetic held. The net energies are computed as follows (ILawson et al.ll2015l) : 


TME 

PME 

NAME 

TKE 

NAKE 


1 fO.l /•7r/2 f2Tx 

III I —2 2 
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1 fTr/2 r^n _ _ 
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2po 

1 

2/io 


'0.61 JO JO 

pO.7 ^7r/2 p27r 

/o .61 Jo J 


y^r"^ sin 9dcf)d9dr , 



( 11 ) 


where TME, PME and TKE, are the axi-symmetric toroidal and poloidal magnetic en¬ 
ergies and the azimuthal kinetic energy, respectively. NAME and NAKE stand for non- 
axisymmetric magnetic and kinetic energies. The resulting evolution prohles are presented 
in Fig. |9l The axi-symmetric toroidal (poloidal) magnetic and kinetic energies are depicted 
by the continuous black (green) and red lines, the non-axisymmetric magnetic and kinetic 
energies are shown by the dashed blue and orange lines. The results indicate that all models 
develop non-axisymmetric motions and magnetic held in the stable layer, however the energy 
distribution is diherent from model to model. Models RCOl and RC02 reach steady state 
after few years of evolution. In the hrst case the toroidal magnetic held is weak and erratic 
in the radiative layer (see right panel of Figs. [6^ and|9^). The weakness of this component of 
the held in the stable layer probably results from the inhibition of penetrative convection by 
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the fast rotation. The poloidal held has larger energy, which is consistent with the presence 
of a non-zero a effect in this region (red line in Fig. [3 ]d). 


The case of model RC02 is more interesting since there is an important mean oscillatory 
toroidal held. The non-axisymmetric magnetic held is also oscillatory and exhibits a phase lag 
with respect to the toroidal held (see the vertical dashed lines indicating particnlar peak time s 
for TME and NAME ). This behavionr is similar to the one obtained bv iMiesch et al.l (120071 ): 


Lawson et al.l (120151) who observe the growth and decay of these qnantities. Nonetheless, in 


these works the valnes of NAME and TME are similar to each other and the phase lag is 
abont 7t/2. In the model RC02 the toroidal energy is larger than the tnrbnlent magnetic 
energy. It is noteworthy the fact that the mean poloidal held oscillates with a 7r/2 phase 
lag with respect to the non-axisymmetric magnetic held (i.e., when the PME grows the 
NAME decays). Note also that the PME is smaller than the NAME which is snggestive of 
a magnetic a-ehect giving rise to a large-scale poloidal magnetic held. The energy of the 
tnrbnlent motions oscillates with the same phase than the tnrbnlent held bnt with lower 
amplitude. Finally, note also that the mean kinetic energy is larger than the toroidal energy, 
the hgure indicates acc eleration and decelera tion of the angular velocity in phase with the 
toroidal magnetic held (iGuerrero et al.ll2016l) . 


The model RC03 is the one that takes longer to reach steady state. Since this case 
is the one that exhibits more variation, in Fig. |9t we present the entire evolution of the 
energies in the stable layer. The hgure shows that the small scale non-axisymmetric kinetic 
and magnetic held energy components develop after the development of the toroidal held and 
that the evolution of these two quantities is closely tight (see dashed lines). The mean toroidal 
kinetic energy evolves together with the mean toroidal held. When the non-axisymmetric 
modes develop, the toroidal kinetic energy grows. We interpret this behaviour as the inwards 
transport of angular momentum, which makes the stable region to accelerate. The poloidal 
magnetic energy also develops after the establishment of MHD turbulence in the radiative 
zone. After around 350 years of evolution the rotation seems to achieve a steady state. This 
modihes all the other quantities, specially the mean poloidal energy, which grows about 2 
orders of magnitude in few years. Note that the a-ehect has important positive values in 
the entire stable layer (see the prohle of a for the model RC03, green line in Fig. 13b which 
was computed in the hnal phase of the evolution). 


In general, it is evident that the dynamics of the convectively stable stratihed layer 
increases in complexity for models with larger Rossby number. This, of course, means a 
dependence on the unstable region where phenomena like penetrative convection or turbulent 
pumping are established. The non-axisymmetric how and held seem to be contributing to 
the generation of poloidal held. From the FOSA analyses presented above, it is the small 
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scale magnetic field which is contributing through the magnetic a-effect. The turbulent flows 
in the radiative zone increase the turbulent diffusion in this region. For the oscillatory model 
RC02, it dehnes the cycle period of the large-scale magnetic held in the entire domain. 



33 


a) RCOl - Energy [10^° J] 



time [years] 

b) RC02 - Energy [10^" J] 




100 200 300 400 

time [years] 


Fig. 9.— Volume integral of the different energies in the stable region for the models a) 
RCOl, b) RC02 and c) RC03). The toroidal magnetic energy (TME) is shown in black, the 
toroidal kinetic energy (TKE) in red, the poloidal held energy (PEE) in green, the non- 
axisymmetric magnetic held (NAME) with the dashed blue line and the non-axysymmetric 
kinetic energy (NAKE) with the dashed orange line. The dotted lines in the middle panels 
highlight the phase-lag between the axi-symmetric and the non-axisymmetric quantities. In 
the bottom panel, the dotted line shows the point from which the bntterhy diagram of Fig. 
Et was plotted. 
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4. Summary and conclusions 


We have performed global 3D MHD simulations of turbulent convection in rotating 
spherical shells using the anelastic approximation. We have used the code EULAG-MHD 
which captures the contribution of the non-resolved scales via an implicit sub-grid scale 
model. The main goal of this work was the comparison of the properties of global dynamos 
for models of the convective envelope only (CZ models), and models that include a part of 
the radiative core below the convection zone (RC models). The stratihcation is set in such a 
way that the Rossby number (Ro) is roughly the same in corresponding CZ and RC models. 
These simulations aimed at contributing to a broader understanding of the amplihcation and 
dynamics of large-scale magnetic helds in the Sun and solar-type stars. Our main conclusions 
are summarized and briefly discussed below. 


The resulting differential rotation prohle depends on the balance between the Coriolis, 
buoyancy and Lorentz forces. All the cases presented here develop a solar-like differential 
rotation (with faster rotation at the equator and slower rotation at the poles). This result 



Nevertheless, our results are in agreement with recent global models that obtain the 
transition from anti-solar to solar-like rotation due to the presence of dynarn o-generated 
magnetic helds flFan &: Fandl2014J: ISimitev et al.ll2015[) . Also. iKarak et al.l (120151) have stud¬ 
ied the dynamo action in a broad range of Rossby numbers and found an anti-solar rota¬ 
tion pattern in slowly rotating stars. In our simulations, the slowest rotating models with 
D = D 0 / 2 , CZ03 and RC03, develop a large-scale magnetic held whose energy is several 
orders of magnitude smaller than the kinetic energy. We did not obtain a regime at which 
the rotation pattern switches from solar to anti-solar. 


The models without a tachocline develop substantial latitudinal diherential rotation. 
Tilted, conical, iso-rotation contours are obtained for the models with relatively slower rota¬ 
tion (i.e., models CZ02 and CZ03). The results reproduce remarkably well the diherential 
rotation in the solar convection zone, with the equator-to-60° latitude diherence, xn ~ 0.18 
(the same as in the Sun), and also reproduce the near-surface shear layer (NSSL). 

Because of the highly sub-adiabatic stable region, the simulations including the upper 
part of the radiative interior develop a tachocline, i.e., a rotational shear layer between the 
radiative and convective zones. The model with faster rotation (D = 2f2Q), RCOl, does not 
develop a NSSL. The slower rotating models, which have a larger inhuence of the buoyancy 
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force, RC02 (f2 = ^q) and RC03 (fi = VLq/2), do develop a NSSL. 

The meridional circulation in the models with radiative zone show a poleward flow at 
latitudes > 30° and a multicellular pattern radially distributed. This last characteristic re¬ 
sembles and theoretically supports recent helioseismology results. It is noteworthy, however, 
that to this date no global simulation has been able to accurately reproduce both, the solar 
differential rotation and the meridional circulation. The computational difficulties to simu¬ 
late the processes at the upper part of the convection zone are perhaps the reason for the 
general lack of success. 


Large-scale dynamo action is observed in all models. The differential rotation, the dy¬ 
namo growth rate and magnetic field topology depend on the Rossby number. The dynamo 
solutions include both steady and oscillatory dynamo regimes. The ratio between the mag¬ 
netic and kinetic energies increases with the rotation rate, i.e. , the smaller Ro the larger 
cm/ck (see Table [T]). Our results partially agree with those of IPetit et al.l (120081) since we 
observe a linear dependence of the mean magnetic energy with the rotation rate. Further¬ 
more, in the fast rotating model with the tachocline, RCOl, the magnetic energy is deposited 
in the poloidal magnetic field. Although most of this poloidal field is concentrated at the 
tachocline, this is reminiscent of the magne tic structure of fast rotating young stars, whose 
observations indicate strong poloidal fields fjCregorv et ah 2012 ). 


Using our simulations results, we have computed the inductive and dissipative terms of 
the mean-field induction equation. We notice first that in both sets of dynamo models, the 
inductive term due to the radial shear is inversely proportional to the rotation rate. The 
shear localizes mainly in the boundary layers, i.e., the interface between the radiative and 
convective zones and in the upper surface. On the other hand, the inductive term due to 
the latitudinal shear in the simulations without the tachocline is proportional to Oq- In fhe 
models RCOl, RC02 and RC03, when the rotation is fast, dgVL peaks at the equator (RCOl). 
When the rotation is slow, dgVL peaks at the poles (RC03). 


The turbulent coefficients and a were computed using the FOSA approximation. 
We observe no dependence of the kinetic counterpart of a, on the Rossby number. 
However, the magnetic counterpart a^n (Eqs. 6 and 7) clearly anti-correlates with Ro; i.e., 
the smaller Ro, the larger the magnetic field and, thus, the larger the small-scale current 
helicity. Consequently, the magnetic a-effect, computed here via the current helicity, could 
work as an inductive term of the magii etic field generation. This partially agrees with the 
results of IVishniac fc Shapovalov! (120141) in the sense that both, kinetic and magnetic parts 
of the a-effect may act as a source of magnetic field. 


In the oscillatory dynamo models without tachocline CZOl and CZ02, the evolution 
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of the field is consistent with dynamo waves propagating according to the Parker-Yoshimura 
sign rule. It is not the case, however, for the tachocline model, RC02, in which the migration 
of the held only in some regions can be explained by this rule. Because the evolution of the 
held in this model is slow (on a time-scale of ~ 30 years), other factors like the meridional 
circulation and/or turbulent magnetic pumping, could be inhuencing the held transport. 


Though the FOSA is the simplest form of computing the mean-held coefficients, in 
general terms our 3D simulation results give support to the mean-held dynamo theory as 
a formalism to explain sol ar and stellar dynamos. A more sophisticate d technique, like 


the test-held method fe.g.. ISchrinner et al.ll2005l: iBrandenburg et al.l 120081) . could provide a 


better inference of these coefficients and thus allow a better analysis. 


The cycle period for the models without tachocline, CZOl and CZ03, is short (about 2 
years) compared to the solar cycle period. In the slow and fast rotating tachocline models, 
RCOl and RC03, the dynamo regimes are stationary. However in the model RC02, with 
the solar rotation rate and the tachocline, the held is oscillatory with the period of about 
30 years. This time scale is comparable with the solar cycle, however, unlike in the Sun, the 
toroidal magnetic held is symmetric across the equator. 


In the RC models, strong toroidal magnetic held generated by dynamo is concentrated 
at or below the interface between the convective and the radiative layers. This indicates 
that the main source of the magnetic held is the radial shear at the tachocline. This ehect 
has enormous inhuence on the general evolution of the models. First, this strong magnetic 
held modihes the mean-how prohles via the Maxwell stresses; and second, because the held 
is deposited in a convectively stable region, its evolution occurs on longer time-scales than 
in the models without the tachocline. Despite the held production is distributed over the 
simulation domain, the dynamo period is regulated by the dihusion time of the toroidal 
magnetic held at the deeper, more stable layers. The deep seated magnetic held seems to 
control the evolution of the magnetic held in the rest of the convection zone. This result 
could explanain a long standing problem in mean-held modeling, i.e., the coexistence of a 
highly turbulent magnetic dihusivity, of the order of 10® m^ s“^, with the 22-year cycle 
period of the Sun. 


The value of the magnetic dihusivity, and so of the dihusion time, in the stable layer is 
a product of the small scale motions generated by MHD instabilities. We have studied the 
temporal evolution of the mean, axi-symmetric, as well as the non-axisymmetric kinetic and 
magnetic energy components. This temporal evolution suggests that once the toroidal held is 
established due to the strong radial shear, magneto-shear instabilities develop in the radiative 
zone generating turbulent motions and magnetic helds. The results support the idea of a 
magnetic a-ehect as a source of poloidal held. The turbulent motions, on the other hand. 
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enhance the magnetic diffusivity. Therefore, the system adjusts itself to the new quantities 
until a new steady state is reached. The dynamics in the stable layer depends on the Rossby 
number, slow rotating models exhibit more complex behaviour and take longer to achieve 
steady state. In model RC02 the turbulent moti ons in the stab l e layer define the ~ 30-year 
cycle period. The solar model of tachocline of iGhizaru et al.l (1201011 uses a stratification 
that results in slower convective, less penetrative, motions. For that reason the cycle period 
is longer that the one found here. The difference in the cycle period between models with 
and without tachocline, as well as the develpment of a magnetic a-effect, suggest that the 
instabilities in the tachocline might be of significant importance in the solar dynamo. 


One could ask why RC models do not exhibit the same short-period dynamo that 
is observed in the CZ models due to the latidudinal shear alone. The answer relies on the 
backreaction of the magnetic field on the fluid dynamics via the Lorentz force. This influence 
can be easily noticed by comparing the profiles of O between the models CZ02 and RC02. In 
numerical experiments, not presented here, we have verified that the strength of the toroidal 
magnetic field (and thus its influence on the flow) at the tachocline can be controlled by 
varying different model ingredients. The radial shear source term is particularly sensitive to 
the ambient state profile, ©e. For instance, changes in parameters like the thickness of the 
transition between radiative and convective zones, the amplitude of the thermal oscillations 
or the adiabaticity of radiative and convective layers, could result in different rotation and 
dynamo patterns. The solar values of these parameters are not fully established for the Sun. 
We should also stress that a better understanding of the effects of the turbulent diffusion. 


de Gouveia Dal Pino et ah 

2012 ; 

Karak et ah 

2014; 

Santos-Lima et ah 

2015) 


We conclude that tachoclines play an important role in solar and stellar dynamos. 
Modelers should be careful when interpreting results obtained in simulations that do not 
include this rotational shear layer. From this work, several new studies seem to be promising 
for the understanding of the solar/stellar rotation and magnetism. For instance, a study 
of mean-flows and angular momentum balance for cool stars in the presence of magnetic 
fields could help to understand the cyclic and non-cyclic stellar activity. Models without 
tachocline, which evolve in short time-scales, and could be interpreted more easily, could 
help in the study of the origin of stochasticity and intermittency in turbulent flows. This 
will provide new insights on the physical origin of Mounder-like solar minima. The results 
are also encouraging to the study of the rotation evolution and magnetism in young stars. 
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